# "Media's Influence on LGBTQ Support Across Africa" by Stephen Winkler

# Code to generate figure with placebo results:
#   - Figure 3 in main text

# Content:
# - Load and transform data 
# - Run  model for each independent variable (radio, tv, newspaper, internet, social media)
#   on each dependent variable (religion, ethnicity, lgbtq, hiv, immigrant)
# - Extract point estimats & vcov in each model 

rm(list=ls()) # clear environment
source("Rprofile.R") # setup Rprofile (load packages, etc.)

# Load data 
myData <- readRDS("data/afrobarometer.rds") #if using rds file

# Convert Country into Dummy & Rename
dummies <- predict(dummyVars(~ country, data = myData, levelsOnly = TRUE), newdata = myData)
colnames(dummies) <- gsub("country", "", colnames(dummies)) #remove country from column names
colnames(dummies) <- gsub(" ", "", colnames(dummies)) #remove spaces
colnames(dummies) <- gsub("'", "", colnames(dummies)) # remove unusual characters

# Join country dummies back to main dataset
myData <- cbind(myData, dummies)
myData <- dplyr::select(myData, -Algeria, -Sudan, -Egypt) #remove NA countries

# Transform data to numeric b/c sims doesn't like logical/factors
myData$radio <- as.numeric(myData$radio)
myData$tv <- as.numeric(myData$tv)
myData$newspaper <- as.numeric(myData$newspaper)
myData$internet <- as.numeric(myData$internet)
myData$social_media <- as.numeric(myData$social_media)
myData$media_aggregate <- as.numeric(myData$media_aggregate)
myData$media_noradio <- as.numeric(myData$media_noradio)
myData$media_notv <- as.numeric(myData$media_notv)
myData$media_nopaper <- as.numeric(myData$media_nopaper)
myData$media_nointernet <- as.numeric(myData$media_nointernet)
myData$media_nosocialmedia <- as.numeric(myData$media_nosocialmedia)
myData$sexuality2 <- as.numeric(myData$sexuality2)
myData$relig_tol2 <- as.numeric(myData$relig_tol2)
myData$ethnic_tol2 <- as.numeric(myData$ethnic_tol2)
myData$hiv_tol2 <- as.numeric(myData$hiv_tol2)
myData$immig_tol2 <- as.numeric(myData$immig_tol2)
myData$tolerance_noLGBT2 <- as.numeric(myData$tolerance_noLGBT2)
myData$female <- as.numeric(myData$female)
myData$education <- as.numeric(myData$education)
myData$religiosity <- as.numeric(myData$religiosity)
myData$age <- as.character(myData$age) #age needs to go through character first
myData$age <- as.numeric(myData$age) #then to numeric to maintain real values
myData$income_water <- as.numeric(myData$income_water)
myData$urban <- as.numeric(myData$urban)

## Media Aggregate

# religion
model.media_aggregate.relig <- 
  relig_tol2 ~
  media_aggregate  +
  ethnic_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.media_aggregate.relig <- extractdata(model.media_aggregate.relig, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
media_aggregate.relig.result <- glm(model.media_aggregate.relig, family = binomial(link="logit"), data = mdata.media_aggregate.relig) #run logit
vc.media_aggregate.relig <- cluster.vcov(media_aggregate.relig.result, mdata.media_aggregate.relig$district) #extract clustered vcov matrix
media_aggregate.relig.result.robust <- coeftest(media_aggregate.relig.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.media_aggregate.relig <- media_aggregate.relig.result.robust[,1] #extract pe from clustered model 

# ethnicity
model.media_aggregate.ethn <- 
  ethnic_tol2 ~
  media_aggregate + 
  relig_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.media_aggregate.ethn <- extractdata(model.media_aggregate.ethn, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
media_aggregate.ethn.result <- glm(model.media_aggregate.ethn, family = binomial(link="logit"), data = mdata.media_aggregate.ethn) #run logit
vc.media_aggregate.ethn <- cluster.vcov(media_aggregate.ethn.result, mdata.media_aggregate.ethn$district) #extract clustered vcov matrix
media_aggregate.ethn.result.robust <- coeftest(media_aggregate.ethn.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.media_aggregate.ethn <- media_aggregate.ethn.result.robust[,1] #extract pe from clustered model 

# lbtq
model.media_aggregate.lgbtq <- 
  sexuality2 ~
  media_aggregate + 
  tolerance_noLGBT2 + #add social tolerance 
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.media_aggregate.lgbtq <- extractdata(model.media_aggregate.lgbtq, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
media_aggregate.lgbtq.result <- glm(model.media_aggregate.lgbtq, family = binomial(link="logit"), data = mdata.media_aggregate.lgbtq) #run logit
vc.media_aggregate.lgbtq <- cluster.vcov(media_aggregate.lgbtq.result, mdata.media_aggregate.lgbtq$district) #extract clustered vcov matrix
media_aggregate.lgbtq.result.robust <- coeftest(media_aggregate.lgbtq.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.media_aggregate.lgbtq <- media_aggregate.lgbtq.result.robust[,1] #extract pe from clustered model 

# hiv
model.media_aggregate.hiv <- 
  hiv_tol2 ~
  media_aggregate + 
  ethnic_tol2 + sexuality2 + relig_tol2 + immig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.media_aggregate.hiv <- extractdata(model.media_aggregate.hiv, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
media_aggregate.hiv.result <- glm(model.media_aggregate.hiv, family = binomial(link="logit"), data = mdata.media_aggregate.hiv) #run logit
vc.media_aggregate.hiv <- cluster.vcov(media_aggregate.hiv.result, mdata.media_aggregate.hiv$district) #extract clustered vcov matrix
media_aggregate.hiv.result.robust <- coeftest(media_aggregate.hiv.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.media_aggregate.hiv <- media_aggregate.hiv.result.robust[,1] #extract pe from clustered model 

# immigrants
model.media_aggregate.immig <- 
  immig_tol2 ~
  media_aggregate +
  ethnic_tol2 + sexuality2 + hiv_tol2 + relig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.media_aggregate.immig <- extractdata(model.media_aggregate.immig, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
media_aggregate.immig.result <- glm(model.media_aggregate.immig, family = binomial(link="logit"), data = mdata.media_aggregate.immig) #run logit
vc.media_aggregate.immig <- cluster.vcov(media_aggregate.immig.result, mdata.media_aggregate.immig$district) #extract clustered vcov matrix
media_aggregate.immig.result.robust <- coeftest(media_aggregate.immig.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.media_aggregate.immig <- media_aggregate.immig.result.robust[,1] #extract pe from clustered model 

## RADIO

# religion
model.radio.relig <- 
  relig_tol2 ~
  radio + media_noradio +
  ethnic_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.radio.relig <- extractdata(model.radio.relig, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
radio.relig.result <- glm(model.radio.relig, family = binomial(link="logit"), data = mdata.radio.relig) #run logit
vc.radio.relig <- cluster.vcov(radio.relig.result, mdata.radio.relig$district) #extract clustered vcov matrix
radio.relig.result.robust <- coeftest(radio.relig.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.radio.relig <- radio.relig.result.robust[,1] #extract pe from clustered model 

# ethnicity
model.radio.ethn <- 
  ethnic_tol2 ~
  radio + media_noradio +
  relig_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.radio.ethn <- extractdata(model.radio.ethn, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
radio.ethn.result <- glm(model.radio.ethn, family = binomial(link="logit"), data = mdata.radio.ethn) #run logit
vc.radio.ethn <- cluster.vcov(radio.ethn.result, mdata.radio.ethn$district) #extract clustered vcov matrix
radio.ethn.result.robust <- coeftest(radio.ethn.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.radio.ethn <- radio.ethn.result.robust[,1] #extract pe from clustered model 

# lbtq
model.radio.lgbtq <- 
  sexuality2 ~
  radio + media_noradio +
  tolerance_noLGBT2 + #add social tolerance 
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.radio.lgbtq <- extractdata(model.radio.lgbtq, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
radio.lgbtq.result <- glm(model.radio.lgbtq, family = binomial(link="logit"), data = mdata.radio.lgbtq) #run logit
vc.radio.lgbtq <- cluster.vcov(radio.lgbtq.result, mdata.radio.lgbtq$district) #extract clustered vcov matrix
radio.lgbtq.result.robust <- coeftest(radio.lgbtq.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.radio.lgbtq <- radio.lgbtq.result.robust[,1] #extract pe from clustered model 

# hiv
model.radio.hiv <- 
  hiv_tol2 ~
  radio + media_noradio +
  ethnic_tol2 + sexuality2 + relig_tol2 + immig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.radio.hiv <- extractdata(model.radio.hiv, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
radio.hiv.result <- glm(model.radio.hiv, family = binomial(link="logit"), data = mdata.radio.hiv) #run logit
vc.radio.hiv <- cluster.vcov(radio.hiv.result, mdata.radio.hiv$district) #extract clustered vcov matrix
radio.hiv.result.robust <- coeftest(radio.hiv.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.radio.hiv <- radio.hiv.result.robust[,1] #extract pe from clustered model 

# immigrants
model.radio.immig <- 
  immig_tol2 ~
  radio + media_noradio +
  ethnic_tol2 + sexuality2 + hiv_tol2 + relig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.radio.immig <- extractdata(model.radio.immig, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
radio.immig.result <- glm(model.radio.immig, family = binomial(link="logit"), data = mdata.radio.immig) #run logit
vc.radio.immig <- cluster.vcov(radio.immig.result, mdata.radio.immig$district) #extract clustered vcov matrix
radio.immig.result.robust <- coeftest(radio.immig.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.radio.immig <- radio.immig.result.robust[,1] #extract pe from clustered model 

## TV
# religion
model.tv.relig <- 
  relig_tol2 ~
  tv + media_notv +
  ethnic_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.tv.relig <- extractdata(model.tv.relig, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
tv.relig.result <- glm(model.tv.relig, family = binomial(link="logit"), data = mdata.tv.relig) #run logit
vc.tv.relig <- cluster.vcov(tv.relig.result, mdata.tv.relig$district) #extract clustered vcov matrix
tv.relig.result.robust <- coeftest(tv.relig.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.tv.relig <- tv.relig.result.robust[,1] #extract pe from clustered model 

# ethnicity
model.tv.ethn <- 
  ethnic_tol2 ~
  tv + media_notv +
  relig_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.tv.ethn <- extractdata(model.tv.ethn, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
tv.ethn.result <- glm(model.tv.ethn, family = binomial(link="logit"), data = mdata.tv.ethn) #run logit
vc.tv.ethn <- cluster.vcov(tv.ethn.result, mdata.tv.ethn$district) #extract clustered vcov matrix
tv.ethn.result.robust <- coeftest(tv.ethn.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.tv.ethn <- tv.ethn.result.robust[,1] #extract pe from clustered model 

# lbtq
model.tv.lgbtq <- 
  sexuality2 ~
  tv + media_notv +
  tolerance_noLGBT2 + #add social tolerance 
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.tv.lgbtq <- extractdata(model.tv.lgbtq, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
tv.lgbtq.result <- glm(model.tv.lgbtq, family = binomial(link="logit"), data = mdata.tv.lgbtq) #run logit
vc.tv.lgbtq <- cluster.vcov(tv.lgbtq.result, mdata.tv.lgbtq$district) #extract clustered vcov matrix
tv.lgbtq.result.robust <- coeftest(tv.lgbtq.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.tv.lgbtq <- tv.lgbtq.result.robust[,1] #extract pe from clustered model 

# hiv
model.tv.hiv <- 
  hiv_tol2 ~
  tv + media_notv +
  ethnic_tol2 + sexuality2 + relig_tol2 + immig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.tv.hiv <- extractdata(model.tv.hiv, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
tv.hiv.result <- glm(model.tv.hiv, family = binomial(link="logit"), data = mdata.tv.hiv) #run logit
vc.tv.hiv <- cluster.vcov(tv.hiv.result, mdata.tv.hiv$district) #extract clustered vcov matrix
tv.hiv.result.robust <- coeftest(tv.hiv.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.tv.hiv <- tv.hiv.result.robust[,1] #extract pe from clustered model 

# immigrants
model.tv.immig <- 
  immig_tol2 ~
  tv + media_notv +
  ethnic_tol2 + sexuality2 + hiv_tol2 + relig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.tv.immig <- extractdata(model.tv.immig, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
tv.immig.result <- glm(model.tv.immig, family = binomial(link="logit"), data = mdata.tv.immig) #run logit
vc.tv.immig <- cluster.vcov(tv.immig.result, mdata.tv.immig$district) #extract clustered vcov matrix
tv.immig.result.robust <- coeftest(tv.immig.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.tv.immig <- tv.immig.result.robust[,1] #extract pe from clustered model 

## NEWSPAPER
# religion
model.newspaper.relig <- 
  relig_tol2 ~
  newspaper + media_nopaper +
  ethnic_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.newspaper.relig <- extractdata(model.newspaper.relig, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
newspaper.relig.result <- glm(model.newspaper.relig, family = binomial(link="logit"), data = mdata.newspaper.relig) #run logit
vc.newspaper.relig <- cluster.vcov(newspaper.relig.result, mdata.newspaper.relig$district) #extract clustered vcov matrix
newspaper.relig.result.robust <- coeftest(newspaper.relig.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.newspaper.relig <- newspaper.relig.result.robust[,1] #extract pe from clustered model 

# ethnicity
model.newspaper.ethn <- 
  ethnic_tol2 ~
  newspaper + media_nopaper +
  relig_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.newspaper.ethn <- extractdata(model.newspaper.ethn, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
newspaper.ethn.result <- glm(model.newspaper.ethn, family = binomial(link="logit"), data = mdata.newspaper.ethn) #run logit
vc.newspaper.ethn <- cluster.vcov(newspaper.ethn.result, mdata.newspaper.ethn$district) #extract clustered vcov matrix
newspaper.ethn.result.robust <- coeftest(newspaper.ethn.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.newspaper.ethn <- newspaper.ethn.result.robust[,1] #extract pe from clustered model 

# lbtq
model.newspaper.lgbtq <- 
  sexuality2 ~
  newspaper + media_nopaper +
  tolerance_noLGBT2 + #add social tolerance 
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.newspaper.lgbtq <- extractdata(model.newspaper.lgbtq, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
newspaper.lgbtq.result <- glm(model.newspaper.lgbtq, family = binomial(link="logit"), data = mdata.newspaper.lgbtq) #run logit
vc.newspaper.lgbtq <- cluster.vcov(newspaper.lgbtq.result, mdata.newspaper.lgbtq$district) #extract clustered vcov matrix
newspaper.lgbtq.result.robust <- coeftest(newspaper.lgbtq.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.newspaper.lgbtq <- newspaper.lgbtq.result.robust[,1] #extract pe from clustered model 

# hiv
model.newspaper.hiv <- 
  hiv_tol2 ~
  newspaper + media_nopaper +
  ethnic_tol2 + sexuality2 + relig_tol2 + immig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.newspaper.hiv <- extractdata(model.newspaper.hiv, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
newspaper.hiv.result <- glm(model.newspaper.hiv, family = binomial(link="logit"), data = mdata.newspaper.hiv) #run logit
vc.newspaper.hiv <- cluster.vcov(newspaper.hiv.result, mdata.newspaper.hiv$district) #extract clustered vcov matrix
newspaper.hiv.result.robust <- coeftest(newspaper.hiv.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.newspaper.hiv <- newspaper.hiv.result.robust[,1] #extract pe from clustered model 

# immigrants
model.newspaper.immig <- 
  immig_tol2 ~
  newspaper + media_nopaper +
  ethnic_tol2 + sexuality2 + hiv_tol2 + relig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.newspaper.immig <- extractdata(model.newspaper.immig, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
newspaper.immig.result <- glm(model.newspaper.immig, family = binomial(link="logit"), data = mdata.newspaper.immig) #run logit
vc.newspaper.immig <- cluster.vcov(newspaper.immig.result, mdata.newspaper.immig$district) #extract clustered vcov matrix
newspaper.immig.result.robust <- coeftest(newspaper.immig.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.newspaper.immig <- newspaper.immig.result.robust[,1] #extract pe from clustered model 

## INTERNET
# religion
model.internet.relig <- 
  relig_tol2 ~
  internet + media_nointernet +
  ethnic_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.internet.relig <- extractdata(model.internet.relig, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
internet.relig.result <- glm(model.internet.relig, family = binomial(link="logit"), data = mdata.internet.relig) #run logit
vc.internet.relig <- cluster.vcov(internet.relig.result, mdata.internet.relig$district) #extract clustered vcov matrix
internet.relig.result.robust <- coeftest(internet.relig.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.internet.relig <- internet.relig.result.robust[,1] #extract pe from clustered model 

# ethnicity
model.internet.ethn <- 
  ethnic_tol2 ~
  internet + media_nointernet +
  relig_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.internet.ethn <- extractdata(model.internet.ethn, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
internet.ethn.result <- glm(model.internet.ethn, family = binomial(link="logit"), data = mdata.internet.ethn) #run logit
vc.internet.ethn <- cluster.vcov(internet.ethn.result, mdata.internet.ethn$district) #extract clustered vcov matrix
internet.ethn.result.robust <- coeftest(internet.ethn.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.internet.ethn <- internet.ethn.result.robust[,1] #extract pe from clustered model 

# lbtq
model.internet.lgbtq <- 
  sexuality2 ~
  internet + media_nointernet +
  tolerance_noLGBT2 + #add social tolerance 
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.internet.lgbtq <- extractdata(model.internet.lgbtq, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
internet.lgbtq.result <- glm(model.internet.lgbtq, family = binomial(link="logit"), data = mdata.internet.lgbtq) #run logit
vc.internet.lgbtq <- cluster.vcov(internet.lgbtq.result, mdata.internet.lgbtq$district) #extract clustered vcov matrix
internet.lgbtq.result.robust <- coeftest(internet.lgbtq.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.internet.lgbtq <- internet.lgbtq.result.robust[,1] #extract pe from clustered model 

# hiv
model.internet.hiv <- 
  hiv_tol2 ~
  internet + media_nointernet +
  ethnic_tol2 + sexuality2 + relig_tol2 + immig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.internet.hiv <- extractdata(model.internet.hiv, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
internet.hiv.result <- glm(model.internet.hiv, family = binomial(link="logit"), data = mdata.internet.hiv) #run logit
vc.internet.hiv <- cluster.vcov(internet.hiv.result, mdata.internet.hiv$district) #extract clustered vcov matrix
internet.hiv.result.robust <- coeftest(internet.hiv.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.internet.hiv <- internet.hiv.result.robust[,1] #extract pe from clustered model 

# immigrants
model.internet.immig <- 
  immig_tol2 ~
  internet + media_nointernet +
  ethnic_tol2 + sexuality2 + hiv_tol2 + relig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.internet.immig <- extractdata(model.internet.immig, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
internet.immig.result <- glm(model.internet.immig, family = binomial(link="logit"), data = mdata.internet.immig) #run logit
vc.internet.immig <- cluster.vcov(internet.immig.result, mdata.internet.immig$district) #extract clustered vcov matrix
internet.immig.result.robust <- coeftest(internet.immig.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.internet.immig <- internet.immig.result.robust[,1] #extract pe from clustered model 

# Social media
# religion
model.social_media.relig <- 
  relig_tol2 ~
  social_media + media_nosocialmedia +
  ethnic_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.social_media.relig <- extractdata(model.social_media.relig, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
social_media.relig.result <- glm(model.social_media.relig, family = binomial(link="logit"), data = mdata.social_media.relig) #run logit
vc.social_media.relig <- cluster.vcov(social_media.relig.result, mdata.social_media.relig$district) #extract clustered vcov matrix
social_media.relig.result.robust <- coeftest(social_media.relig.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.social_media.relig <- social_media.relig.result.robust[,1] #extract pe from clustered model 

# ethnicity
model.social_media.ethn <- 
  ethnic_tol2 ~
  social_media + media_nosocialmedia +
  relig_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.social_media.ethn <- extractdata(model.social_media.ethn, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
social_media.ethn.result <- glm(model.social_media.ethn, family = binomial(link="logit"), data = mdata.social_media.ethn) #run logit
vc.social_media.ethn <- cluster.vcov(social_media.ethn.result, mdata.social_media.ethn$district) #extract clustered vcov matrix
social_media.ethn.result.robust <- coeftest(social_media.ethn.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.social_media.ethn <- social_media.ethn.result.robust[,1] #extract pe from clustered model 

# lbtq
model.social_media.lgbtq <- 
  sexuality2 ~
  social_media + media_nosocialmedia +
  tolerance_noLGBT2 + #add social tolerance 
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.social_media.lgbtq <- extractdata(model.social_media.lgbtq, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
social_media.lgbtq.result <- glm(model.social_media.lgbtq, family = binomial(link="logit"), data = mdata.social_media.lgbtq) #run logit
vc.social_media.lgbtq <- cluster.vcov(social_media.lgbtq.result, mdata.social_media.lgbtq$district) #extract clustered vcov matrix
social_media.lgbtq.result.robust <- coeftest(social_media.lgbtq.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.social_media.lgbtq <- social_media.lgbtq.result.robust[,1] #extract pe from clustered model 

# hiv
model.social_media.hiv <- 
  hiv_tol2 ~
  social_media + media_nosocialmedia +
  ethnic_tol2 + sexuality2 + relig_tol2 + immig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.social_media.hiv <- extractdata(model.social_media.hiv, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
social_media.hiv.result <- glm(model.social_media.hiv, family = binomial(link="logit"), data = mdata.social_media.hiv) #run logit
vc.social_media.hiv <- cluster.vcov(social_media.hiv.result, mdata.social_media.hiv$district) #extract clustered vcov matrix
social_media.hiv.result.robust <- coeftest(social_media.hiv.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.social_media.hiv <- social_media.hiv.result.robust[,1] #extract pe from clustered model 

# immigrants
model.social_media.immig <- 
  immig_tol2 ~
  social_media + media_nosocialmedia +
  ethnic_tol2 + sexuality2 + hiv_tol2 + relig_tol2 + # social tolerance of other groups
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata.social_media.immig <- extractdata(model.social_media.immig, myData, na.rm = TRUE, extra = ~respno + district) #only extract data used in model
social_media.immig.result <- glm(model.social_media.immig, family = binomial(link="logit"), data = mdata.social_media.immig) #run logit
vc.social_media.immig <- cluster.vcov(social_media.immig.result, mdata.social_media.immig$district) #extract clustered vcov matrix
social_media.immig.result.robust <- coeftest(social_media.immig.result, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
pe.social_media.immig <- social_media.immig.result.robust[,1] #extract pe from clustered model 

# Initialize scenarios to mean values of covariates
media_aggregate.scen.relig <- cfMake(model.media_aggregate.relig, data=mdata.media_aggregate.relig, nscen=1) #media_aggregate
media_aggregate.scen.ethn <- cfMake(model.media_aggregate.ethn, data=mdata.media_aggregate.ethn, nscen=1)
media_aggregate.scen.lgbtq <- cfMake(model.media_aggregate.lgbtq, data=mdata.media_aggregate.lgbtq, nscen=1)
media_aggregate.scen.hiv <- cfMake(model.media_aggregate.hiv, data=mdata.media_aggregate.hiv, nscen=1)
media_aggregate.scen.immig <- cfMake(model.media_aggregate.immig, data=mdata.media_aggregate.immig, nscen=1)
radio.scen.relig <- cfMake(model.radio.relig, data=mdata.radio.relig, nscen=1) #RADIO
radio.scen.ethn <- cfMake(model.radio.ethn, data=mdata.radio.ethn, nscen=1)
radio.scen.lgbtq <- cfMake(model.radio.lgbtq, data=mdata.radio.lgbtq, nscen=1)
radio.scen.hiv <- cfMake(model.radio.hiv, data=mdata.radio.hiv, nscen=1)
radio.scen.immig <- cfMake(model.radio.immig, data=mdata.radio.immig, nscen=1)
tv.scen.relig <- cfMake(model.tv.relig, data=mdata.tv.relig, nscen=1)           #TV
tv.scen.ethn <- cfMake(model.tv.ethn, data=mdata.tv.ethn, nscen=1)
tv.scen.lgbtq <- cfMake(model.tv.lgbtq, data=mdata.tv.lgbtq, nscen=1)
tv.scen.hiv <- cfMake(model.tv.hiv, data=mdata.tv.hiv, nscen=1)
tv.scen.immig <- cfMake(model.tv.immig, data=mdata.tv.immig, nscen=1)
newspaper.scen.relig <- cfMake(model.newspaper.relig, data=mdata.newspaper.relig, nscen=1) #NEWSPAPER
newspaper.scen.ethn <- cfMake(model.newspaper.ethn, data=mdata.newspaper.ethn, nscen=1)
newspaper.scen.lgbtq <- cfMake(model.newspaper.lgbtq, data=mdata.newspaper.lgbtq, nscen=1)
newspaper.scen.hiv <- cfMake(model.newspaper.hiv, data=mdata.newspaper.hiv, nscen=1)
newspaper.scen.immig <- cfMake(model.newspaper.immig, data=mdata.newspaper.immig, nscen=1)
internet.scen.relig <- cfMake(model.internet.relig, data=mdata.internet.relig, nscen=1)   #INTERNET
internet.scen.ethn <- cfMake(model.internet.ethn, data=mdata.internet.ethn, nscen=1)
internet.scen.lgbtq <- cfMake(model.internet.lgbtq, data=mdata.internet.lgbtq, nscen=1)
internet.scen.hiv <- cfMake(model.internet.hiv, data=mdata.internet.hiv, nscen=1)
internet.scen.immig <- cfMake(model.internet.immig, data=mdata.internet.immig, nscen=1)
social_media.scen.relig <- cfMake(model.social_media.relig, data=mdata.social_media.relig, nscen=1) #SOCIAL MEDIA
social_media.scen.ethn <- cfMake(model.social_media.ethn, data=mdata.social_media.ethn, nscen=1)
social_media.scen.lgbtq <- cfMake(model.social_media.lgbtq, data=mdata.social_media.lgbtq, nscen=1)
social_media.scen.hiv <- cfMake(model.social_media.hiv, data=mdata.social_media.hiv, nscen=1)
social_media.scen.immig <- cfMake(model.social_media.immig, data=mdata.social_media.immig, nscen=1)

# Configure scenarios:  Media Aggregate from 1 to 5 
media_aggregate.scen.relig <- cfName(media_aggregate.scen.relig, "media_aggregate 1 vs 21", scen=1) # name the change
media_aggregate.scen.relig <- cfChange(media_aggregate.scen.relig, "media_aggregate", x = 21, xpre=1, scen=1) # set the pre and post
media_aggregate.scen.ethn <- cfName(media_aggregate.scen.ethn, "media_aggregate 1 vs 21", scen=1)
media_aggregate.scen.ethn <- cfChange(media_aggregate.scen.ethn, "media_aggregate", x = 21, xpre=1, scen=1)
media_aggregate.scen.lgbtq <- cfName(media_aggregate.scen.lgbtq, "media_aggregate 1 vs 21", scen=1)
media_aggregate.scen.lgbtq <- cfChange(media_aggregate.scen.lgbtq, "media_aggregate", x = 21, xpre=1, scen=1)
media_aggregate.scen.hiv <- cfName(media_aggregate.scen.hiv, "media_aggregate 1 vs 21", scen=1)
media_aggregate.scen.hiv <- cfChange(media_aggregate.scen.hiv, "media_aggregate", x = 21, xpre=1, scen=1)
media_aggregate.scen.immig <- cfName(media_aggregate.scen.immig, "media_aggregate 1 vs 21", scen=1)
media_aggregate.scen.immig <- cfChange(media_aggregate.scen.immig, "media_aggregate", x = 21, xpre=1, scen=1)

# Configure scenarios:  RADIO from 1 to 5 
radio.scen.relig <- cfName(radio.scen.relig, "radio 1 vs 5", scen=1) # name the change
radio.scen.relig <- cfChange(radio.scen.relig, "radio", x = 5, xpre=1, scen=1) # set the pre and post
radio.scen.ethn <- cfName(radio.scen.ethn, "radio 1 vs 5", scen=1)
radio.scen.ethn <- cfChange(radio.scen.ethn, "radio", x = 5, xpre=1, scen=1)
radio.scen.lgbtq <- cfName(radio.scen.lgbtq, "radio 1 vs 5", scen=1)
radio.scen.lgbtq <- cfChange(radio.scen.lgbtq, "radio", x = 5, xpre=1, scen=1)
radio.scen.hiv <- cfName(radio.scen.hiv, "radio 1 vs 5", scen=1)
radio.scen.hiv <- cfChange(radio.scen.hiv, "radio", x = 5, xpre=1, scen=1)
radio.scen.immig <- cfName(radio.scen.immig, "radio 1 vs 5", scen=1)
radio.scen.immig <- cfChange(radio.scen.immig, "radio", x = 5, xpre=1, scen=1)

# Configure scenarios:  TV from 1 to 5 
tv.scen.relig <- cfName(tv.scen.relig, "tv 1 vs 5", scen=1) # name the change
tv.scen.relig <- cfChange(tv.scen.relig, "tv", x = 5, xpre=1, scen=1) # set the pre and post
tv.scen.ethn <- cfName(tv.scen.ethn, "tv 1 vs 5", scen=1)
tv.scen.ethn <- cfChange(tv.scen.ethn, "tv", x = 5, xpre=1, scen=1)
tv.scen.lgbtq <- cfName(tv.scen.lgbtq, "tv 1 vs 5", scen=1)
tv.scen.lgbtq <- cfChange(tv.scen.lgbtq, "tv", x = 5, xpre=1, scen=1)
tv.scen.hiv <- cfName(tv.scen.hiv, "tv 1 vs 5", scen=1)
tv.scen.hiv <- cfChange(tv.scen.hiv, "tv", x = 5, xpre=1, scen=1)
tv.scen.immig <- cfName(tv.scen.immig, "tv 1 vs 5", scen=1)
tv.scen.immig <- cfChange(tv.scen.immig, "tv", x = 5, xpre=1, scen=1)

# Configure scenarios:  NEWSPAPER from 1 to 5 
newspaper.scen.relig <- cfName(newspaper.scen.relig, "newspaper 1 vs 5", scen=1) # name the change
newspaper.scen.relig <- cfChange(newspaper.scen.relig, "newspaper", x = 5, xpre=1, scen=1) # set the pre and post
newspaper.scen.ethn <- cfName(newspaper.scen.ethn, "newspaper 1 vs 5", scen=1)
newspaper.scen.ethn <- cfChange(newspaper.scen.ethn, "newspaper", x = 5, xpre=1, scen=1)
newspaper.scen.lgbtq <- cfName(newspaper.scen.lgbtq, "newspaper 1 vs 5", scen=1)
newspaper.scen.lgbtq <- cfChange(newspaper.scen.lgbtq, "newspaper", x = 5, xpre=1, scen=1)
newspaper.scen.hiv <- cfName(newspaper.scen.hiv, "newspaper 1 vs 5", scen=1)
newspaper.scen.hiv <- cfChange(newspaper.scen.hiv, "newspaper", x = 5, xpre=1, scen=1)
newspaper.scen.immig <- cfName(newspaper.scen.immig, "newspaper 1 vs 5", scen=1)
newspaper.scen.immig <- cfChange(newspaper.scen.immig, "newspaper", x = 5, xpre=1, scen=1)

# Configure scenarios:  INTERNET from 1 to 5 
internet.scen.relig <- cfName(internet.scen.relig, "internet 1 vs 5", scen=1) # name the change
internet.scen.relig <- cfChange(internet.scen.relig, "internet", x = 5, xpre=1, scen=1) # set the pre and post
internet.scen.ethn <- cfName(internet.scen.ethn, "internet 1 vs 5", scen=1)
internet.scen.ethn <- cfChange(internet.scen.ethn, "internet", x = 5, xpre=1, scen=1)
internet.scen.lgbtq <- cfName(internet.scen.lgbtq, "internet 1 vs 5", scen=1)
internet.scen.lgbtq <- cfChange(internet.scen.lgbtq, "internet", x = 5, xpre=1, scen=1)
internet.scen.hiv <- cfName(internet.scen.hiv, "internet 1 vs 5", scen=1)
internet.scen.hiv <- cfChange(internet.scen.hiv, "internet", x = 5, xpre=1, scen=1)
internet.scen.immig <- cfName(internet.scen.immig, "internet 1 vs 5", scen=1)
internet.scen.immig <- cfChange(internet.scen.immig, "internet", x = 5, xpre=1, scen=1)

# Configure scenarios:  SOCIAL MEDIA from 1 to 5 
social_media.scen.relig <- cfName(social_media.scen.relig, "social_media 1 vs 5", scen=1) # name the change
social_media.scen.relig <- cfChange(social_media.scen.relig, "social_media", x = 5, xpre=1, scen=1) # set the pre and post
social_media.scen.ethn <- cfName(social_media.scen.ethn, "social_media 1 vs 5", scen=1)
social_media.scen.ethn <- cfChange(social_media.scen.ethn, "social_media", x = 5, xpre=1, scen=1)
social_media.scen.lgbtq <- cfName(social_media.scen.lgbtq, "social_media 1 vs 5", scen=1)
social_media.scen.lgbtq <- cfChange(social_media.scen.lgbtq, "social_media", x = 5, xpre=1, scen=1)
social_media.scen.hiv <- cfName(social_media.scen.hiv, "social_media 1 vs 5", scen=1)
social_media.scen.hiv <- cfChange(social_media.scen.hiv, "social_media", x = 5, xpre=1, scen=1)
social_media.scen.immig <- cfName(social_media.scen.immig, "social_media 1 vs 5", scen=1)
social_media.scen.immig <- cfChange(social_media.scen.immig, "social_media", x = 5, xpre=1, scen=1)

# Simulate conditional expectations for these counterfactuals
sims <- 10000

# Run simulations and generate quantities of interest
# Media Aggregate
simbetas.media_aggregate.relig <- mvrnorm(sims, pe.media_aggregate.relig, vc.media_aggregate.relig) #religion
media_aggregate.relig.qoi <- logitsimfd(media_aggregate.scen.relig, simbetas.media_aggregate.relig, ci=0.95)
simbetas.media_aggregate.ethn <- mvrnorm(sims, pe.media_aggregate.ethn, vc.media_aggregate.ethn) #ethnicity
media_aggregate.ethn.qoi <- logitsimfd(media_aggregate.scen.ethn, simbetas.media_aggregate.ethn, ci=0.95)
simbetas.media_aggregate.lgbtq <- mvrnorm(sims, pe.media_aggregate.lgbtq, vc.media_aggregate.lgbtq) #lgbtq
media_aggregate.lgbtq.qoi <- logitsimfd(media_aggregate.scen.lgbtq, simbetas.media_aggregate.lgbtq, ci=0.95)
simbetas.media_aggregate.hiv <- mvrnorm(sims, pe.media_aggregate.hiv, vc.media_aggregate.hiv) #hiv
media_aggregate.hiv.qoi <- logitsimfd(media_aggregate.scen.hiv, simbetas.media_aggregate.hiv, ci=0.95)
simbetas.media_aggregate.immig <- mvrnorm(sims, pe.media_aggregate.immig, vc.media_aggregate.immig) #immigrants
media_aggregate.immig.qoi <- logitsimfd(media_aggregate.scen.immig, simbetas.media_aggregate.immig, ci=0.95)
# Radio
simbetas.radio.relig <- mvrnorm(sims, pe.radio.relig, vc.radio.relig) #religion
radio.relig.qoi <- logitsimfd(radio.scen.relig, simbetas.radio.relig, ci=0.95)
simbetas.radio.ethn <- mvrnorm(sims, pe.radio.ethn, vc.radio.ethn) #ethnicity
radio.ethn.qoi <- logitsimfd(radio.scen.ethn, simbetas.radio.ethn, ci=0.95)
simbetas.radio.lgbtq <- mvrnorm(sims, pe.radio.lgbtq, vc.radio.lgbtq) #lgbtq
radio.lgbtq.qoi <- logitsimfd(radio.scen.lgbtq, simbetas.radio.lgbtq, ci=0.95)
simbetas.radio.hiv <- mvrnorm(sims, pe.radio.hiv, vc.radio.hiv) #hiv
radio.hiv.qoi <- logitsimfd(radio.scen.hiv, simbetas.radio.hiv, ci=0.95)
simbetas.radio.immig <- mvrnorm(sims, pe.radio.immig, vc.radio.immig) #immigrants
radio.immig.qoi <- logitsimfd(radio.scen.immig, simbetas.radio.immig, ci=0.95)
# TV
simbetas.tv.relig <- mvrnorm(sims, pe.tv.relig, vc.tv.relig) #religion
tv.relig.qoi <- logitsimfd(tv.scen.relig, simbetas.tv.relig, ci=0.95)
simbetas.tv.ethn <- mvrnorm(sims, pe.tv.ethn, vc.tv.ethn) #ethnicity
tv.ethn.qoi <- logitsimfd(tv.scen.ethn, simbetas.tv.ethn, ci=0.95)
simbetas.tv.lgbtq <- mvrnorm(sims, pe.tv.lgbtq, vc.tv.lgbtq) #lgbtq
tv.lgbtq.qoi <- logitsimfd(tv.scen.lgbtq, simbetas.tv.lgbtq, ci=0.95)
simbetas.tv.hiv <- mvrnorm(sims, pe.tv.hiv, vc.tv.hiv) #hiv
tv.hiv.qoi <- logitsimfd(tv.scen.hiv, simbetas.tv.hiv, ci=0.95)
simbetas.tv.immig <- mvrnorm(sims, pe.tv.immig, vc.tv.immig) #immigrants
tv.immig.qoi <- logitsimfd(tv.scen.immig, simbetas.tv.immig, ci=0.95)
# Newspaper
simbetas.newspaper.relig <- mvrnorm(sims, pe.newspaper.relig, vc.newspaper.relig) #religion
newspaper.relig.qoi <- logitsimfd(newspaper.scen.relig, simbetas.newspaper.relig, ci=0.95)
simbetas.newspaper.ethn <- mvrnorm(sims, pe.newspaper.ethn, vc.newspaper.ethn) #ethnicity
newspaper.ethn.qoi <- logitsimfd(newspaper.scen.ethn, simbetas.newspaper.ethn, ci=0.95)
simbetas.newspaper.lgbtq <- mvrnorm(sims, pe.newspaper.lgbtq, vc.newspaper.lgbtq) #lgbtq
newspaper.lgbtq.qoi <- logitsimfd(newspaper.scen.lgbtq, simbetas.newspaper.lgbtq, ci=0.95)
simbetas.newspaper.hiv <- mvrnorm(sims, pe.newspaper.hiv, vc.newspaper.hiv) #hiv
newspaper.hiv.qoi <- logitsimfd(newspaper.scen.hiv, simbetas.newspaper.hiv, ci=0.95)
simbetas.newspaper.immig <- mvrnorm(sims, pe.newspaper.immig, vc.newspaper.immig) #immigrants
newspaper.immig.qoi <- logitsimfd(newspaper.scen.immig, simbetas.newspaper.immig, ci=0.95)
# Internet
simbetas.internet.relig <- mvrnorm(sims, pe.internet.relig, vc.internet.relig) #religion
internet.relig.qoi <- logitsimfd(internet.scen.relig, simbetas.internet.relig, ci=0.95)
simbetas.internet.ethn <- mvrnorm(sims, pe.internet.ethn, vc.internet.ethn) #ethnicity
internet.ethn.qoi <- logitsimfd(internet.scen.ethn, simbetas.internet.ethn, ci=0.95)
simbetas.internet.lgbtq <- mvrnorm(sims, pe.internet.lgbtq, vc.internet.lgbtq) #lgbtq
internet.lgbtq.qoi <- logitsimfd(internet.scen.lgbtq, simbetas.internet.lgbtq, ci=0.95)
simbetas.internet.hiv <- mvrnorm(sims, pe.internet.hiv, vc.internet.hiv) #hiv
internet.hiv.qoi <- logitsimfd(internet.scen.hiv, simbetas.internet.hiv, ci=0.95)
simbetas.internet.immig <- mvrnorm(sims, pe.internet.immig, vc.internet.immig) #immigrants
internet.immig.qoi <- logitsimfd(internet.scen.immig, simbetas.internet.immig, ci=0.95)
# Social Media
simbetas.social_media.relig <- mvrnorm(sims, pe.social_media.relig, vc.social_media.relig) #religion
social_media.relig.qoi <- logitsimfd(social_media.scen.relig, simbetas.social_media.relig, ci=0.95)
simbetas.social_media.ethn <- mvrnorm(sims, pe.social_media.ethn, vc.social_media.ethn) #ethnicity
social_media.ethn.qoi <- logitsimfd(social_media.scen.ethn, simbetas.social_media.ethn, ci=0.95)
simbetas.social_media.lgbtq <- mvrnorm(sims, pe.social_media.lgbtq, vc.social_media.lgbtq) #lgbtq
social_media.lgbtq.qoi <- logitsimfd(social_media.scen.lgbtq, simbetas.social_media.lgbtq, ci=0.95)
simbetas.social_media.hiv <- mvrnorm(sims, pe.social_media.hiv, vc.social_media.hiv) #hiv
social_media.hiv.qoi <- logitsimfd(social_media.scen.hiv, simbetas.social_media.hiv, ci=0.95)
simbetas.social_media.immig <- mvrnorm(sims, pe.social_media.immig, vc.social_media.immig) #immigrants
social_media.immig.qoi <- logitsimfd(social_media.scen.immig, simbetas.social_media.immig, ci=0.95)

# Convert to dataframe so I can plot with ggplot
# media_aggregate
media_aggregate.relig.df <- as.data.frame(media_aggregate.relig.qoi) %>% mutate(medium = "Media aggregate", dv = "Other Religion", model = "logit")
media_aggregate.ethn.df <- as.data.frame(media_aggregate.ethn.qoi) %>% mutate(medium = "Media aggregate", dv = "Other Ethnicity", model = "logit")
media_aggregate.lgbtq.df <- as.data.frame(media_aggregate.lgbtq.qoi) %>% mutate(medium = "Media aggregate", dv = "Homosexuality", model = "logit")
media_aggregate.hiv.df <- as.data.frame(media_aggregate.hiv.qoi) %>% mutate(medium = "Media aggregate", dv = "HIV+", model = "logit")
media_aggregate.immig.df <- as.data.frame(media_aggregate.immig.qoi) %>% mutate(medium = "Media aggregate", dv = "Immigrant", model = "logit")
# radio
radio.relig.df <- as.data.frame(radio.relig.qoi) %>% mutate(medium = "Radio", dv = "Other Religion", model = "logit")
radio.ethn.df <- as.data.frame(radio.ethn.qoi) %>% mutate(medium = "Radio", dv = "Other Ethnicity", model = "logit")
radio.lgbtq.df <- as.data.frame(radio.lgbtq.qoi) %>% mutate(medium = "Radio", dv = "Homosexuality", model = "logit")
radio.hiv.df <- as.data.frame(radio.hiv.qoi) %>% mutate(medium = "Radio", dv = "HIV+", model = "logit")
radio.immig.df <- as.data.frame(radio.immig.qoi) %>% mutate(medium = "Radio", dv = "Immigrant", model = "logit")
# tv
tv.relig.df <- as.data.frame(tv.relig.qoi) %>% mutate(medium = "TV", dv = "Other Religion", model = "logit")
tv.ethn.df <- as.data.frame(tv.ethn.qoi) %>% mutate(medium = "TV", dv = "Other Ethnicity", model = "logit")
tv.lgbtq.df <- as.data.frame(tv.lgbtq.qoi) %>% mutate(medium = "TV", dv = "Homosexuality", model = "logit")
tv.hiv.df <- as.data.frame(tv.hiv.qoi) %>% mutate(medium = "TV", dv = "HIV+", model = "logit")
tv.immig.df <- as.data.frame(tv.immig.qoi) %>% mutate(medium = "TV", dv = "Immigrant", model = "logit")
# newspaper
newspaper.relig.df <- as.data.frame(newspaper.relig.qoi) %>% mutate(medium = "Newspaper", dv = "Other Religion", model = "logit")
newspaper.ethn.df <- as.data.frame(newspaper.ethn.qoi) %>% mutate(medium = "Newspaper", dv = "Other Ethnicity", model = "logit")
newspaper.lgbtq.df <- as.data.frame(newspaper.lgbtq.qoi) %>% mutate(medium = "Newspaper", dv = "Homosexuality", model = "logit")
newspaper.hiv.df <- as.data.frame(newspaper.hiv.qoi) %>% mutate(medium = "Newspaper", dv = "HIV+", model = "logit")
newspaper.immig.df <- as.data.frame(newspaper.immig.qoi) %>% mutate(medium = "Newspaper", dv = "Immigrant", model = "logit")
# internet
internet.relig.df <- as.data.frame(internet.relig.qoi) %>% mutate(medium = "Internet", dv = "Other Religion", model = "logit")
internet.ethn.df <- as.data.frame(internet.ethn.qoi) %>% mutate(medium = "Internet", dv = "Other Ethnicity", model = "logit")
internet.lgbtq.df <- as.data.frame(internet.lgbtq.qoi) %>% mutate(medium = "Internet", dv = "Homosexuality", model = "logit")
internet.hiv.df <- as.data.frame(internet.hiv.qoi) %>% mutate(medium = "Internet", dv = "HIV+", model = "logit")
internet.immig.df <- as.data.frame(internet.immig.qoi) %>% mutate(medium = "Internet", dv = "Immigrant", model = "logit")
# social_media
social_media.relig.df <- as.data.frame(social_media.relig.qoi) %>% mutate(medium = "Social media", dv = "Other Religion", model = "logit")
social_media.ethn.df <- as.data.frame(social_media.ethn.qoi) %>% mutate(medium = "Social media", dv = "Other Ethnicity", model = "logit")
social_media.lgbtq.df <- as.data.frame(social_media.lgbtq.qoi) %>% mutate(medium = "Social media", dv = "Homosexuality", model = "logit")
social_media.hiv.df <- as.data.frame(social_media.hiv.qoi) %>% mutate(medium = "Social media", dv = "HIV+", model = "logit")
social_media.immig.df <- as.data.frame(social_media.immig.qoi) %>% mutate(medium = "Social media", dv = "Immigrant", model = "logit")

all_qoi_robust_placebo <- bind_rows(media_aggregate.relig.df, media_aggregate.ethn.df, media_aggregate.lgbtq.df, media_aggregate.hiv.df, media_aggregate.immig.df,
                                    radio.relig.df, radio.ethn.df, radio.lgbtq.df, radio.hiv.df, radio.immig.df,
                                    tv.relig.df, tv.ethn.df, tv.lgbtq.df, tv.hiv.df, tv.immig.df,
                                    newspaper.relig.df, newspaper.ethn.df, newspaper.lgbtq.df, newspaper.hiv.df, newspaper.immig.df,
                                    internet.relig.df, internet.ethn.df, internet.lgbtq.df, internet.hiv.df, internet.immig.df,
                                    social_media.relig.df, social_media.ethn.df, social_media.lgbtq.df, social_media.hiv.df, social_media.immig.df)

# PLOT
all_qoi_robust_placebo$dv <- factor(all_qoi_robust_placebo$dv,  # set as factor so I can create levels
                                    levels = c("Homosexuality",
                                               "HIV+",
                                               "Other Religion",
                                               "Other Ethnicity",
                                               "Immigrant")) 
all_qoi_robust_placebo$medium <- factor(all_qoi_robust_placebo$medium,
                                        levels = c("Radio",
                                                   "TV",
                                                   "Newspaper",
                                                   "Internet",
                                                   "Social media", 
                                                   "Media aggregate")) 
all_qoi_robust_placebo <- all_qoi_robust_placebo %>% group_by(medium) #group so the qoi of each medium don't get split up
dodge <- position_dodge(width =.4) #set offset positioning

### FIGURE 3 in Main Text:
placebo_plot_logit <- ggplot(all_qoi_robust_placebo) + 
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper, 
                      x = (medium),
                      color = (dv),
                      group = (dv),
                      shape = (dv)),
                  position = dodge) + 
  labs(x = "", y = "difference in probability of supporting out-group (95% CI)") +
  #ggtitle("Change in out-group support when moving from \n 'none' to 'daily' media consumption (logit)") + 
  coord_flip() + 
  #scale_x_discrete(labels=c("radio", "tV", "newspaper", "internet", "social_media")) + 
  scale_y_continuous(labels = scales::percent_format(accuracy = 2), breaks = c(-0.03, -0.02, -0.01, 0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, .07, .08, .09, .1)) +
  scale_color_manual(name = "Out Group",
                     values = c("black", "gray70", "gray70", "gray70", "gray70")) + 
                     #values = rev(brewer.pal(5, "Paired"))) +
  scale_shape_manual(name = "Out Group",
                     values = rev(c(0,1,2,3,8))) +
  #ylim(-0.06,0.06, labels=scales::percent) +
  geom_hline(yintercept=c(0), linetype="dotted") + 
  theme_bw() + theme(panel.border=element_blank(), 
                     panel.grid.minor = element_blank(),
                     panel.grid.major = element_blank(),
                     plot.title = element_text(hjust=0.5),
                     legend.position = c(.9,.5))
placebo_plot_logit

ggsave(filename = "figures/placebo_plot_logit.pdf", #Save it to directory
       plot = placebo_plot_logit,
       width = 8, height = 6.5)